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LONG-TERM  GOALS 

The  long-term  goal  is  to  develop  a  nonhydrostatic,  parallel  ocean  simulation  tool  that  is  capable  of 
simulating  processes  on  a  wide  range  of  scales  through  use  of  accurate  numerical  methods  and  high- 
performance  computational  algorithms.  The  tool  will  be  applied  to  study  highly  nonlinear  internal 
waves  in  coastal  domains  to  develop  an  improved  understanding  of  mechanisms  that  govern  their 
generation,  propagation,  and  dissipation. 

OBJECTIVES 

The  primary  objective  is  to  enhance  the  capabilities  of  the  SUNTANS  model  through  development  of 
algorithms  to  study  multiscale  processes  in  estuaries  and  the  coastal  ocean.  This  involves  development 
of  1)  improved  momentum  and  scalar  advection  on  unstructured,  staggered  grids,  2)  accurate  and 
efficient  algorithms  for  solution  of  the  nonhydrostatic  pressure,  and  3)  adaptive  grid  capabiliites  with 
adaptive  mesh  refinement  and  model  nesting. 

APPROACH 

This  work  focuses  on  the  continued  development  of  SUNTANS  (Stanford  Unstructured 
Nonhydrostatic  Terrain-following  Adaptive  Navier-Stokes  Simulator),  a  free-surface,  nonhydrostatic, 
unstructured-grid,  parallel  coastal  ocean  and  estuary  simulation  tool  that  solves  the  Navier-Stokes 
equations  under  the  Boussinesq  approximation  (Fringer  et  al.,2006).  The  formulation  is  based  on  the 
method  outlined  by  Casulli  and  Walters  (2000),  in  which  the  free-surface  and  vertical  diffusion  are 
discretized  with  the  theta  method  which  eliminates  the  Courant  condition  associated  with  fast  free- 
surface  waves  and  the  elevated  friction  term  associated  with  small  vertical  grid  spacing  at  the  free- 
surface  and  bottom  boundary.  For  flows  with  extensive  wetting  and  drying,  advection  of  momentum  is 
accomplished  with  the  semi-Lagrangian  advection  scheme  (Wang  et  al.  2011a),  which  ensures  stability 
in  the  presence  of  cells  that  fill  and  empty  with  the  tides.  Scalar  advection  is  accomplished  semi- 
implicitly  and  continuity  of  volume  and  mass  are  guaranteed  for  the  hydrostatic  solver.  The  theta 
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method  for  the  free  surface  yields  a  two-dimensional  Poisson  equation,  and  the  nonhydrostatic 
pressure  is  governed  by  a  three-dimensional  Poisson  equation.  These  are  both  solved  with  the 
preconditioned  conjugate  gradient  algorithm  with  Jacobi  and  block-Jacobi  preconditioning, 
respectively.  Because  the  nonhydrostatic  component  of  SUNTANS  is  essentially  a  correction  to  the 
hydrostatic  component,  SUNTANS  can  be  run  seamlessly  in  nonhydrostatic  or  hydrostatic  modes. 
SUNTANS  is  written  in  the  C  programming  language,  and  the  message-passing  interface  (MPI)  is 
employed  for  use  in  a  distributed-memory  parallel  computing  environment.  SUNTANS  employs  the 
generalized  length  scale  approach  to  Reynolds-averaged  turbulence  modeling  (Wang  et  al.  2011b). 

The  SUNTANS  grid  employs  z-levels  in  the  vertical  and  is  unstructured  in  plan,  which  enables  the 
resolution  of  complex  coastlines  and  topographic  features.  Unstructured  grids  also  enable  the  use  of 
high  grid  resolution  in  regions  of  interest  while  coarsening  the  grid  in  regions  where  grid  resolution  is 
not  required,  thereby  significantly  reducing  computational  overhead. 

WORK  COMPLETED 

We  have  completed  the  development  of  a  nonhydrostatic  isopycnal-coordinate  model  and  have 
performed  simulations  to  study  the  nonlinear  effects  of  internal  wave  generation  over  an  idealized 
Gaussian  ridge.  We  have  also  performed  high-resolution  simulations  of  breaking  internal  waves  on 
slopes  and  computed  mixing  and  dissipation.  Below  results  are  presented  for  the  isopycnal-coordinate 
model  and  the  internal  wave  generation. 

RESULTS 

Nonhydrostatic  isopycnal-coordinate  modeling 

We  have  nearly  completed  a  manuscript  reporting  the  formulation  and  testing  of  a  nonhydrostatic 
isopycnal  coordinate  ocean  model,  which  to  our  knowledge  represents  the  first  of  its  kind.  The 
motivation  for  this  model  comes  from  the  success  and  efficiency  of  traditional  (hydrostatic)  isopycnal- 
coordinate  models.  Isopycnal  coordinates  naturally  represent  a  stratified  fluid,  which  in  turn  reduces 
the  number  of  grid  points  from  0(100)  grid  points  in  traditional  ocean  models  (as  in  Zhang  et  al.  2011) 
to  0(1)  grid  points.  This  represents  a  significant  reduction  in  computational  effort.  Traditional 
(hydrostatic)  isopycnal-coordinate  models  have  been  applied  to  simulating  internal  waves  in  the  South 
China  Sea  (Simmons  et  al.  2011).  However,  hydrostatic  models  don’t  compute  physical  dispersion 
and  thus  modeled  solitary  waves  are  non-physical  (they  arise  from  a  balance  between  nonlinearity  and 
spurious  numerical  dispersion,  as  discussed  in  Vitousek  and  Fringer  2011).  We  have  developed  a 
nonhydrostatic  isopycnal-coordinate  ocean  model  to  overcome  this  limitation  and  to  allow  the 
formation  of  solitary  waves  with  an  isopycnal  coordinate  formulation.  Ideally,  in  locations  where  the 
internal  wave  structure  is  predominantly  mode-1,  an  isopycnal  model  with  just  two  layers  may  suffice. 

Our  nonhydrostatic  isopycnal-coordinate  model  is  intended  to  be  flexible  (using  an  arbitrary  number  of 
layers)  and  straightforward  (resembling  existing  ocean  models).  We  employ  a  second-order  accurate 
(in  space  and  time)  finite-difference/finite-volume  formulation  on  a  staggered  C-grid  as  in  traditional 
isopycnal-coordinate  models.  However,  our  model  uses  an  implicit  free-surface  discretization  based 
on  the  IMEX  multistep  methods  of  Durran  and  Blossey  (2012)  rather  than  traditional  mode-splitting. 
We  note  that  implicit  free-surface  discretizations  are  more  commonly  used  in  nonhydrostatic  ocean 
models.  Our  model  uses  the  MPDATA  algorithm  (Smolarkiewicz  and  Margolin  1998)  to  handle 
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wetting  and  drying  of  isopycnal  layers.  The  model  solves  the  nonhydrostatic  pressure  Poisson 
equation  using  either  the  preconditioned  conjugate  gradient  method  or  the  multigrid  algorithm  with 
semi-coarsening  in  the  horizontal  and  line  relaxation  in  the  vertical,  a  strategy  that  performs  well  for 
solving  anisotropic  elliptic  equations  (Briggs  et  al.  2000). 


Idealized  test  cases  show  that  the  model  is  capable  of  resolving  dispersive  wave  properties  and  the 
formation  of  nonlinear  internal  solitary  waves.  Figure  1  shows  some  of  the  arbitrary  layer 
configurations  that  can  be  used  in  the  model,  while  Figure  2  shows  the  dispersive  properties  of  the 
hydrostatic  and  nonhydrostatic  isopycnal  models  which  are  validated  against  linear  theory.  These 
results  suggest  that  nonhydrostatic  effects  require  more  layers  even  if  the  stratification  is  uniform  in 
order  to  resolve  depth-variability  associated  with  nonhydrostatic  effects.  Figure  3  shows  the 
comparison  of  the  formation  of  nonlinear  internal  solitary  waves  in  the  SUNTANS  (z-level  model; 
Fringer  et  al.  2006)  to  the  nonhydrostatic  isopycnal  model  with  2  and  10  layers.  This  test  case 
demonstrates  that  the  nonhydrostatic  isopycnal  model  is  capable  of  representing  the  formation  and 
propagation  of  solitary  waves.  Finally,  we  demonstrate  that,  although  the  model  represents  a  discrete 
set  of  layers,  it  is  capable  of  simulating  internal  wave  behavior  in  continuously  stratified  systems. 
Figure  4  shows  a  nonhydrostatic  isopycnal  simulation  of  oscillatory  flow  in  a  continuously  stratified 
fluid  over  a  Gaussian  sill.  This  depth  perturbation  induces  the  generation  of  internal  wave  beams 
which  radiate  away  from  the  sill  at  an  angle  that  agrees  well  with  the  linear  theory. 


1  layer  4  layers  1 0  layers 


X  |m}  X  [iti]  X  [m] 

Figure  1:  The  arbitrary  layer  configuration  used  in  the  nonhydrostatic  isopycnal  models. 
Panels  A,B  and  C  depict  a  free  surface  seiche,  panels  D,E  and  F  depict  a  2-layer  internal  seiche 
and  panels  G,H  and  I  depict  an  internal  seiche  with  a  hyperbolic  tangent  density  profile  with 

varying  number  of  layers. 
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Figure  2:  Modeled  wave  speed  normalized  by  the  true  wave  speed  for  the  hydrostatic  (red)  and 
nonhydrostatic  (blue)  models  as  a  function  of  aspect  ratio  (kH)  for  the  free  surface,  2-layer  and 
smoothly  stratified  internal  seiche  shown  in  panels  A,B  and  C,  respectively.  These  correspond  to 

columns  A,  B,  and  C  in  Figure  1. 
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Figure  3:  Comparison  of  the  formation  of  nonlinear  internal  solitary  waves  in  the  SUNTANS  z- 
level  model  (Fringer  et  al.  2006)  to  our  nonhydrostatic  isopycnal  model  with  2  and  10  layers.  This 
test  case  is  based  on  a  numerical  experiment  in  Vitousek  and  Fringer  (2011). 
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Figure  4:  The  velocity  perturbation  (u  '=u-U)  calculated  by  the  nonhydrostatic  isopycnal 
model  as  induced  by  an  oscillatory  flow  in  a  constantly  stratified  fluid  over  a  Gaussian  sill. 
Internal  wave  beams  radiate  away  from  the  sill  at  an  angle  that  agrees  with  linear  theory. 

We  note  that  the  magnitude  of  the  velocity  perturbation  (and  the  forcing  velocity)  is  kept 

small  so  that  linear  theory  is  valid. 

Nonlinear  and  nonhydrostatic  internal  tide  generation 

The  problem  of  internal  tide  generation  has  conventionally  been  defined  by  three  dimensionless 
parameters,  namely  the  tidal  excursion  parameter  (s2  =  tidal  excursion/  horizontal  topographic  scale), 
the  criticality  parameter  (si  =  topographic  slope/Intemal  wave  beam  slope),  and  the  nondimensional 
obstacle  height  (6  =  topographic  height/total  water  depth).  These  parameters  are  sufficient  to  describe 
the  problem  in  the  linear  limit  {^2 «  1)^  and  analytical  solutions  under  this  approximation  require 
additional  limitations  of  either  Si  «  1  or  5  «  1  or  both.  Recent  advancements  in  analytical  models 
using  Green’s  functions  have  eliminated  these  additional  limitations.  For  example,  the  iTides  model 
of  Echeverri  and  Peacock  (2011)  does  not  require  the  assumptions  of  small  Si  or  5,  although  a  small 
excursion  parameter  is  still  required.  Because  previous  studies  have  largely  been  based  on  linear 
theory  with  the  aforementioned  limitations,  nonlinear  and  nonhydrostatic  effects  on  radiated  energy 
flux  are  largely  unexplored. 

We  carry  out  two-dimensional,  high-resolution  simulations  using  the  SUNTANS  model  to  study 
nonlinear  and  nonhydrostatic  effects  for  realistic  oceanic  flows  with  S2~0.03-0.3.  The  topography 
consists  of  an  idealized  Gaussian  ridge  centered  at  x=0  of  the  form 


fe(x:)  =  #i(jexp(— x^/2cr^), 
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where  ho  is  the  ridge  height  and  o  is  the  width  of  the  Gaussian.  The  domain  is  forced  with  a  barotropic 
pressure  gradient  oscillating  at  the  M2  tidal  frequency  (00  =  1.41x10'^  rad/s).  The  total  water  depth  is 
H=  1000  m  and  the  initial  background  buoyancy  frequency  is  uniform  at  A^=  0.005  rad/s.  We  employ 
a  quadrilateral  grid  (SUNTANS  can  employ  quadrilaterals  or  triangles  or  both),  and  the  finest 
horizontal  resolution  is  employed  over  the  ridge  with  Ax=  200  m  and  linear  stretching  is  employed 
such  that  the  resolution  near  the  boundaries  is  ranges  from  1400  m  to  6000  m  (where  6000  m 
corresponds  to  the  subcritical  slope  with  Si=0.5  and  1400  m  is  the  highest  resolution  at  the  boundaries 
and  corresponds  to  supercritical  slopes  with  Si>  1.75).  A  total  of  50  grid  points  is  used  to  discretize  the 
vertical  direction  and  3800  to  12000  points  are  employed  in  the  horizontal  (the  highest  number  of  grid 
points  are  used  for  the  case  with  si=0.5),  the  time  step  size  is  At=10  s,  and  the  total  simulation  time  is 
10  tidal  periods.  We  investigate  the  effect  of  the  criticality  parameter  varying  from  subcritical  to  very 
supercritical  slopes  and  the  nondimensional  height  parameter  on  the  radiated  energy  flux  for  values  of 
the  excursion  parameter  ranging  from  linear  to  nonlinear  flows.  The  tidal  excursion  parameter  is 
defined  by  S2=uo/(co  o),  where  uq  is  the  barotropic  current  in  deep  water  away  from  the  ridge,  and  we 
vary  uq  but  keep  co  and  a  fixed.  The  ridge  height  is  chosen  to  be  ho=750  m  which  implies  a  tall  ridge 
with  5  =  0.75  and  is  representative  of  areas  with  intense  internal  tide  generation  like  the  Luzon  Strait  in 
the  South  China  Sea  (e.g.  Buijsman  et  al.  2011)  and  Kaena  ridge  at  Hawaii  (e.g.  Carter  et  ah,  2008). 

We  analyze  the  variation  of  the  linear  baroclinic  energy  flux  (i.e.  depth-integrated  and  time-averaged 
p  ’u  ’)  as  a  function  of  the  criticality  parameter.  In  the  figures,  the  theoretical  flux  is  obtained  with  the 
iTides  model  (Echevery  and  Peacock,  2010),  and  results  from  the  linearized  and  hydrostatic 
SUNTANS  model  are  also  computed  for  comparison.  All  results  are  normalized  by  the  theoretical 
maximum  energy  flux  which  is  given  by  the  knife-edge  result  of  St.  Laurent  et  al.  (2003)  that  occurs  in 
the  limit  £2  00  for  finite  S.  To  obtain  linear,  hydrostatic  SUNTANS  results,  the  nonhydrostatic 

pressure  in  SUNTANS  is  eliminated  and  the  nonlinear  terms  in  the  momentum  equation  are  ignored. 

In  the  density  equation,  the  nonlinear  terms  are  eliminated  to  yield  the  linear  density  equation 


dp 

dt 


Ml 

g 


w. 


where  w  is  the  vertical  velocity,  po  is  the  reference  density,  and  g  is  the  gravitational  acceleration. 

As  shown  in  Figure  5,  the  iTides  and  linear  SUNTANS  results  agree  for  all  values  of  S/  and  £2, 
implying  that,  as  expected,  both  models  produce  the  same  results  independent  of  However, 
increasing  the  excursion  parameter  (right-to-left  tiles  in  Figure  5)  produces  marked  deviation  between 
the  nonlinear  and  nonhydrostatic  SUNTANS  result  and  iTides,  particularly  in  the  near-critical  range  81 
=  0.75-1.75.  Figure  6  shows  the  ratio  of  the  energy  flux  computed  by  the  nonlinear  and  nonhydrostatic 
SUNTANS  model  to  the  linear  result  computed  by  iTides,  and  indicates  that  there  is  a  minimum  in  the 
radiated  energy  flux  when  the  slope  is  near-critical.  When  the  slope  is  exactly  critical,  nonlinear 
effects  induce  a  drop  of  85%  in  the  radiated  energy  flux. 
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Figure  5:  Comparison  of  three  different  models  and  how  they  compute  the  effect  of  the  criticality 
parameter  Sj  on  the  baroclinic  energy  flux  radiating  from  a  Gaussian  ridge  for  three  different 
excursion  parameters  (left:  £2=0.03,  middle:  £2=0.1;  right:  £2=0.3)  when  d=0. 75.  The  energy  flux  is 
nondimensionalized  by  the  knife-edge  energy  flux  from  St.  Laurent  et  al.  (2003). 


Figure  6:  Ratio  of  nonlinear  energy  flux  radiating  from  a  Gaussian  ridge  to  the  linear  energy  flux 
computed  by  iTides  as  a  function  of  the  criticality  parameter  £j  when  £2=0.3  and  S=0. 75. 

Figure  7  shows  that  the  drop  in  energy  flux  for  large  exeursion  parameters  and  near-critieal  slopes 
depends  to  great  extent  on  the  height  of  the  ridge.  Both  the  nonlinear  and  nonhydrostatic  SUNTANS 
model  and  the  iTides  model  predict  a  drop  in  the  energy  flux  with  increasing  ridge  height  (with  £]  and 
£2  fixed).  Because  the  iTides  model  employs  a  nonlinear  kinematic  boundary  condition  at  the  bottom 
boundary,  it  is  valid  for  flows  with  arbitrary  £7  and  S,  as  long  as  £2«1.  Therefore,  the  drop  in  energy 
flux  for  large  ridge  heights  predicted  by  the  iTides  model  is  related  to  nonlinear  and  inviscid 
topographic  boundary  effects  and  not  nonlinear  effects  within  the  flow  itself  or  viscous  boundary 
effects  such  as  those  that  give  rise  to  the  nonlinear  energy  flux  terms  (e.g.  Kang  and  Fringer  2012)  or 
to  breaking  and  dissipation.  These  flow-induced  nonlinear  effects  give  rise  to  further  losses  in  the 
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energy  flux  as  indicated  by  the  lower  energy  flux  computed  by  the  SUNTANS  model  as  compared  to 
iTides. 

The  relative  effects  of  nonlinearity  within  the  flow  to  those  computed  by  iTides  can  be  understood 
from  Figure  8,  which  shows  a  plot  of  the  ratio  of  the  SUNTANS-computed  energy  flux  to  that 
computed  by  iTides.  As  the  ridge  height  increases,  the  relative  impact  of  nonlinearity  within  the  flow 
on  the  energy  flux  increases.  When  S=0. 75,  the  flux  computed  by  SUNTANS  is  only  25%  of  that 
computed  by  iTides.  Because  the  iTides  model  produces  a  decreased  energy  flux  for  large  S,  it  is  clear 
that  much  of  the  reduction  in  energy  flux  with  increased  S  is  not  a  result  of  an  increase  in  the  nonlinear 
energy  flux  terms  of  Kang  and  Fringer  (2012).  We  hypothesize  that  the  reduction  is  due  to  redirection 
of  barotropic  energy  to  generation  of  boundary  currents  for  taller  ridges  which  results  in  a  lower 
conversion  of  barotropic-to-baroclinic  energy  flux.  Following  the  work  of  Winters  and  Armi  (2013), 
we  are  studying  the  problem  in  terms  of  newly  defined  inner  flow  variables  which  more  accurately 
represent  the  strongly  nonlinear  flow  over  the  ridge,  rather  than  the  outer  flow  variables  Si,  S2,  and  5 
that  were  used  in  this  study  as  well  as  throughout  the  literature  on  internal  wave  generation. 


Figure  7:  Effect  of  the  nondimensional  ridge  height  S  on  the  energy  flux  radiated  from  a  Gaussian 
Ridge  with  ei=l  and  62=0.3  as  computed  by  the  linear  iTides  model  and  the  SUNTANS  model.  The 
energy  flux  is  normalized  by  the  maximum  knife-edge  value  from  St.  Laurent  et  al.  (2003). 
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Figure  8:  Effect  of  the  nondimensional  ridge  height  Son  the  ratio  of  the  energy  flux  radiated 
from  a  Gaussian  Ridge  as  computed  by  SUNTANS  to  that  computed  by  the  Hides  model 

with  Si=l  and  62=0.3. 


IMPACT/APPLICATIONS 

High-resolution  simulations  using  nonhydrostatic  models  like  SUNTANS  are  crucial  for  understanding 

multiscale  processes  that  are  unresolved,  and  hence  parameterized,  in  larger-scale  ocean  models. 
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